module GOLD_thickness_diffuse
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, July 1999                                      *
!*                                                                     *
!*    The subroutine in this file implements horizontal interface      *
!*  depth diffusion.  In a Z-coordinate model, this would be described *
!*  as Gent-McWilliams diffusion, with an identical interpretation     *
!*  of the diffusion coefficient.  The coefficient is locally limited  *
!*  to guarantee numerical stability, and are vertically uniform.      *
!*                                                                     *
!*  Macros written all in capital letters are defined in GOLD_memory.h.*
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v, vh, vav                               *
!*    j    x ^ x ^ x   At >:  u, uh, uav                               *
!*    j    > o > o >   At o:  h                                        *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use GOLD_checksums, only : hchksum, uchksum, vchksum
use GOLD_diag_mediator, only : post_data, query_averaging_enabled, diag_ptrs
use GOLD_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
use GOLD_error_handler, only : GOLD_error, FATAL, WARNING
use GOLD_EOS, only : calculate_density, calculate_density_derivs
use GOLD_file_parser, only : get_param, log_version, param_file_type
use GOLD_grid, only : ocean_grid_type
use GOLD_interface_heights, only : find_eta
use GOLD_lateral_mixing_coeffs, only : VarMix_CS
use GOLD_variables, only : thermo_var_ptrs
use GOLD_EOS, only : calculate_density
use GOLD_MEKE_types, only : MEKE_type

implicit none ; private

public thickness_diffuse, thickness_diffuse_init, thickness_diffuse_end
public vert_fill_TS

#include <GOLD_memory.h>

type, public :: thickness_diffuse_CS ; private
  real    :: Khth           ! The background interface depth diffusivity in m2 s-1.
  real    :: Khth_Slope_Cff ! Slope dependence coefficient of Khth in m2 s-1.
  real    :: Khth_Min       ! Minimum value of Khth in m2 s-1.
  real    :: Khth_Max       ! Maximum value of Khth in m2 s-1, or 0 for no max.
  real    :: slope_max      ! Slopes steeper than this are limited in some way.
  real    :: kappa_smooth   ! A diffusivity that is used to interpolate more
                            ! sensible values of T & S into thin layers.
  logical :: thickness_diffuse ! If true, interfaces heights are diffused
                             ! with a coefficient of Khth.
  logical :: detangle_interfaces ! If true, add 3-d structured interface height
                            ! diffusivities to horizonally smooth jagged layers.
  real    :: detangle_time  ! If detangle_interfaces is true, this is the
                            ! timescale over which maximally jagged grid-scale
                            ! thickness variations are suppressed.  This must be
                            ! longer than DT, or 0 (the default) to use DT.
  logical :: bulkmixedlayer ! If true, a refined bulk mixed layer is used.
  integer :: nkml           ! The number of layers within the mixed layer.
  logical :: debug          ! If true, write verbose checksums for debugging purposes.
  type(diag_ptrs), pointer :: diag ! A pointer to a structure of shareable
                             ! ocean diagnostic fields.
  real, pointer :: GMwork(:,:) => NULL()  ! Work by thick. diff. in W m-2.
  integer :: id_uhGM = -1, id_vhGM = -1, id_GMwork = -1, id_KH_u = -1, id_KH_v = -1
  integer :: id_KH_u1 = -1, id_KH_v1 = -1
 ! integer :: id_sfn_slope_x = -1, id_sfn_slope_y = -1, id_sfn_x = -1, id_sfn_y = -1
end type thickness_diffuse_CS

contains

subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, MEKE, VarMix, CS)
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(inout) :: h
  real, dimension(NXMEMQ_,NYMEM_,NZ_), intent(inout) :: uhtr
  real, dimension(NXMEM_,NYMEMQ_,NZ_), intent(inout) :: vhtr
  type(thermo_var_ptrs),               intent(in)    :: tv
  real,                                intent(in)    :: dt
  type(ocean_grid_type),               intent(in)    :: G
  type(MEKE_type),                     pointer       :: MEKE
  type(VarMix_CS),                     pointer       :: VarMix
  type(thickness_diffuse_CS),          pointer       :: CS
!    This subroutine does interface depth diffusion.  The fluxes are
!  limited to give positive definiteness, and the diffusivities are
!  limited to guarantee stability.

! Arguments: h - Layer thickness, in m.
!  (in/out)  uhtr - Accumulated zonal mass fluxes in m3.
!  (in/out)  vhtr - Accumulated meridional mass fluxes in m3.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      dt - Time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in)      VarMix - A structure containing a number of fields related to
!                     variable lateral mixing.
!  (in)      MEKE - A structure containing information about the Mesoscale Eddy
!                   Kinetic Energy parameterization; this might be unassociated.
!  (in)      CS - The control structure returned by a previous call to
!                 thickness_diffuse_init.


  real :: e(SZI_(G), SZJ_(G), SZK_(G)+1) ! The heights of the interfaces, relative
                                    ! to mean sea level, in m, pos. upward.
  real :: uhD(SZIQ_(G), SZJ_(G), SZK_(G)) ! uhD & vhD are the diffusive u*h &
  real :: vhD(SZI_(G), SZJQ_(G), SZK_(G)) ! v*h fluxes, in m3 s-1.

  real, dimension(SZIQ_(G), SZJ_(G), SZK_(G)+1) :: &
    KH_u, &       ! The interface height diffusivities in u-columns, in m2 s-1.
    int_slope_u   ! A nondimensional ratio from 0 to 1 that gives the relative
                  ! weighting of the interface slopes to that calculated also
                  ! using density gradients at u points.  The physically correct
                  ! slopes occur at 0, while 1 is used for numerical closures.
  real, dimension(SZI_(G), SZJQ_(G), SZK_(G)+1) :: &
    KH_v, &       ! The interface height diffusivities in v-columns, in m2 s-1.
    int_slope_v   ! A nondimensional ratio from 0 to 1 that gives the relative
                  ! weighting of the interface slopes to that calculated also
                  ! using density gradients at v points.  The physically correct
                  ! slopes occur at 0, while 1 is used for numerical closures.
  real :: KH_u_CFL(SZIQ_(G), SZJ_(G)) ! The maximum stable interface height
  real :: KH_v_CFL(SZI_(G), SZJQ_(G)) ! diffusivities at u & v grid points,
                                      ! in m2 s-1.
  real :: Khth_Loc      ! Locally calculated thickness mixing coefficient (m2/s)
  logical :: use_VarMix, Resoln_scaled
  integer :: i, j, k, is, ie, js, je, nz


  if (.not. ASSOCIATED(CS)) call GOLD_error(FATAL, "GOLD_thickness_diffuse:"// &
         "Module must be initialized before it is used.")
  if ((.not.CS%thickness_diffuse) .or. &
       .not.( CS%Khth > 0.0 .or. associated(VarMix) .or. associated(MEKE) ) ) return

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

  if (ASSOCIATED(MEKE)) then
    if (ASSOCIATED(MEKE%GM_src)) then
      do j=js,je ; do i=is,ie ; MEKE%GM_src(i,j) = 0. ; enddo ; enddo
    endif
  endif

  use_VarMix = .false. ; Resoln_scaled = .false.
  if (Associated(VarMix)) then
    use_VarMix = VarMix%use_variable_mixing
    Resoln_scaled = VarMix%Resoln_scaled_KhTh
  endif

  call find_eta(h, tv, G%g_Earth, G, e, halo_size=1)

  ! Set the diffusivities.
  do j=js,je ; do I=is-1,ie
    Khth_Loc = CS%Khth
    if (use_VarMix) &
      Khth_Loc = Khth_Loc + CS%KHTH_Slope_Cff*VarMix%L2u(I,j)*VarMix%SN_u(I,j)
    if (ASSOCIATED(MEKE)) then
      if (associated(MEKE%Kh)) &
        Khth_Loc = Khth_Loc + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i+1,j))
    endif
    if (Resoln_scaled) &
      Khth_Loc = Khth_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i+1,j))
    if (CS%Khth_Max > 0) then
      Khth_Loc = max(CS%Khth_min, min(Khth_Loc,CS%Khth_Max))
    else
      Khth_Loc = max(CS%Khth_min, Khth_Loc)
    endif
    KH_u_CFL(I,j) = 0.2/(dt*(G%IDXu(I,j)*G%IDXu(I,j) + G%IDYu(I,j)*G%IDYu(I,j)))
    KH_u(I,j,1) = min(KH_u_CFL(I,j), Khth_Loc)
  enddo ; enddo
  do K=2,nz+1 ; do j=js,je ; do I=is-1,ie
    KH_u(I,j,K) = KH_u(I,j,1)
  enddo ; enddo ; enddo

  do J=js-1,je ; do i=is,ie
    Khth_Loc = CS%Khth
    if (use_VarMix) &
      Khth_Loc = Khth_Loc + CS%KHTH_Slope_Cff*VarMix%L2v(i,J)*VarMix%SN_v(i,J)
    if (ASSOCIATED(MEKE)) then
      if (associated(MEKE%Kh)) &
        Khth_Loc = Khth_Loc + MEKE%KhTh_fac*sqrt(MEKE%Kh(i,j)*MEKE%Kh(i,j+1))
    endif
    if (Resoln_scaled) &
      Khth_Loc = Khth_Loc * 0.5*(VarMix%Res_fn_h(i,j) + VarMix%Res_fn_h(i,j+1))
    if (CS%Khth_Max > 0) then
      Khth_Loc = max(CS%Khth_min, min(Khth_Loc,CS%Khth_Max))
    else
      Khth_Loc = max(CS%Khth_min, Khth_Loc)
    endif
    KH_v_CFL(i,J) = 0.2/(dt*(G%IDXv(i,J)*G%IDXv(i,J) + G%IDYv(i,J)*G%IDYv(i,J)))
    KH_v(i,J,1) = min(KH_v_CFL(i,J), Khth_Loc)
  enddo ; enddo
  do K=2,nz+1 ; do J=js-1,je ; do i=is,ie
    KH_v(i,J,K) = KH_v(i,J,1)
  enddo ; enddo ; enddo

  do K=1,nz+1 ; do j=js,je ; do I=is-1,ie ; int_slope_u(I,j,K) = 0.0 ; enddo ; enddo ; enddo
  do K=1,nz+1 ; do J=js-1,je ; do i=is,ie ; int_slope_v(i,J,K) = 0.0 ; enddo ; enddo ; enddo

  if (CS%detangle_interfaces) then
    call add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, CS, &
                           int_slope_u, int_slope_v)
  endif

  if (CS%debug) then
    call uchksum(Kh_u(:,:,:),"Kh_u",G,haloshift=0)
    call vchksum(Kh_v(:,:,:),"Kh_v",G,haloshift=0)
    call uchksum(int_slope_u(:,:,:),"int_slope_u",G,haloshift=0)
    call vchksum(int_slope_v(:,:,:),"int_slope_v",G,haloshift=0)
    call hchksum(h(:,:,:),"thickness_diffuse_1 h",G,haloshift=0)
    call hchksum(e(:,:,:),"thickness_diffuse_1 e",G,haloshift=0)
  endif

  call thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, CS, &
                              int_slope_u, int_slope_v)

  do k=1,nz
    do j=js,je ; do I=is-1,ie
      uhtr(I,j,k) = uhtr(I,j,k) + uhD(I,j,k)*dt
      if (ASSOCIATED(CS%diag%uhGM)) CS%diag%uhGM(I,j,k) = uhD(I,j,k)
    enddo ; enddo
    do J=js-1,je ; do i=is,ie
      vhtr(i,J,k) = vhtr(i,J,k) + vhD(i,J,k)*dt
      if (ASSOCIATED(CS%diag%vhGM)) CS%diag%vhGM(i,J,k) = vhD(i,J,k)
    enddo ; enddo
    do j=js,je ; do i=is,ie
      h(i,j,k) = h(i,j,k) - dt * G%IDXDYh(i,j) * &
          ((uhD(I,j,k) - uhD(I-1,j,k)) + (vhD(i,J,k) - vhD(i,J-1,k)))
      if (h(i,j,k) < G%Angstrom) h(i,j,k) = G%Angstrom
    enddo ; enddo
  enddo

  if (ASSOCIATED(MEKE) .and. ASSOCIATED(VarMix)) then
    if (ASSOCIATED(MEKE%Rd_dx_h) .and. ASSOCIATED(VarMix%Rd_dx_h)) then
      do j=js,je ; do i=is,ie
        MEKE%Rd_dx_h(i,j) = VarMix%Rd_dx_h(i,j)
      enddo ; enddo
    endif
  endif

  if (CS%debug) then
    call uchksum(uhD(:,:,:),"thickness_diffuse uhD",G,haloshift=0)
    call vchksum(vhD(:,:,:),"thickness_diffuse vhD",G,haloshift=0)
    call uchksum(uhtr(:,:,:),"thickness_diffuse uhtr",G,haloshift=0)
    call vchksum(vhtr(:,:,:),"thickness_diffuse vhtr",G,haloshift=0)
    call hchksum(h(:,:,:),"thickness_diffuse h",G,haloshift=0)
  endif

! Offer diagnostic fields for averaging.
  if (query_averaging_enabled(CS%diag)) then
    if (CS%id_uhGM > 0) call post_data(CS%id_uhGM, CS%diag%uhGM, CS%diag)
    if (CS%id_vhGM > 0) call post_data(CS%id_vhGM, CS%diag%vhGM, CS%diag)
    if (CS%id_GMwork > 0) call post_data(CS%id_GMwork, CS%GMwork, CS%diag)
    if (CS%id_KH_u > 0) call post_data(CS%id_KH_u, KH_u, CS%diag)
    if (CS%id_KH_v > 0) call post_data(CS%id_KH_v, KH_v, CS%diag)
    if (CS%id_KH_u1 > 0) call post_data(CS%id_KH_u1, KH_u(:,:,1), CS%diag)
    if (CS%id_KH_v1 > 0) call post_data(CS%id_KH_v1, KH_v(:,:,1), CS%diag)
  endif

end subroutine thickness_diffuse

subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, dt, G, MEKE, &
                                  CS, int_slope_u, int_slope_v)
  real, dimension(NXMEM_,NYMEM_,NZ_),    intent(in)  :: h
  real, dimension(NXMEM_,NYMEM_,NZp1_),  intent(in)  :: e
  real, dimension(NXMEMQ_,NYMEM_,NZp1_), intent(in)  :: Kh_u
  real, dimension(NXMEM_,NYMEMQ_,NZp1_), intent(in)  :: Kh_v
  type(thermo_var_ptrs),                 intent(in)  :: tv
  real, dimension(NXMEMQ_,NYMEM_,NZ_),   intent(out) :: uhD
  real, dimension(NXMEM_,NYMEMQ_,NZ_),   intent(out) :: vhD
  real,                                  intent(in)  :: dt
  type(ocean_grid_type),                 intent(in)  :: G
  type(MEKE_type),                       pointer     :: MEKE
  type(thickness_diffuse_CS),            pointer     :: CS
  real, dimension(NXMEMQ_,NYMEM_,NZp1_), optional, intent(in)  :: int_slope_u
  real, dimension(NXMEM_,NYMEMQ_,NZp1_), optional, intent(in)  :: int_slope_v
!    This subroutine does interface depth diffusion.  The fluxes are
!  limited to give positive definiteness, and the diffusivities are
!  limited to guarantee stability.

! Arguments: h - Layer thickness, in m.
!  (in)      e - Interface heights relative to mean sea level, in m.
!  (in)      Kh_u - Thickness diffusivity on interfaces at u points, in m2 s-1.
!  (in)      Kh_v - Thickness diffusivity on interfaces at v points, in m2 s-1.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (out)     uhD - Zonal mass fluxes in m3 s-1.
!  (out)     vhD - Meridional mass fluxes in m3 s-1.
!  (in)      dt - Time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in)      MEKE - A structure containing information about the Mesoscale Eddy
!                   Kinetic Energy parameterization; this might be unassociated.
!  (in)      CS - The control structure returned by a previous call to
!                 thickness_diffuse_init.
!  (in,opt)  int_slope_u - Ratios that determine how much of the isopycnal
!                 slopes are taken directly from the interface slopes without
!                 consideration of density gradients.
!  (in,opt)  int_slope_v - Ratios that determine how much of the isopycnal
!                 slopes are taken directly from the interface slopes without
!                 consideration of density gradients.

  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
    T, &          ! The temperature (or density) in C, with the values in
                  ! in massless layers filled vertically by diffusion.
    S, &          ! The filled salinity, in PSU, with the values in
                  ! in massless layers filled vertically by diffusion.
    Rho, &        ! Density itself, when a nonlinear equation of state is
                  ! not in use.
    h_avail, &    ! The mass available for diffusion out of each face, divided
                  ! by dt, in m3 s-1.
    h_frac        ! The fraction of the mass in the column above the bottom
                  ! interface of a layer that is within a layer, ND. 0<h_frac<=1
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
    pres, &       ! The pressure at an interface, in Pa.
    h_avail_rsum  ! The running sum of h_avail above an interface, in m3 s-1.
  real, dimension(SZIQ_(G), SZJ_(G)) :: &
    drho_dT_u, &  ! The derivatives of density with temperature and
    drho_dS_u     ! salinity at u points, in kg m-3 K-1 and kg m-3 psu-1.
  real, dimension(SZI_(G), SZJQ_(G)) :: &
    drho_dT_v, &  ! The derivatives of density with temperature and
    drho_dS_v     ! salinity at v points, in kg m-3 K-1 and kg m-3 psu-1.
  real :: uhtot(SZIQ_(G), SZJ_(G))  ! The vertical sum of uhD, in m3 s-1.
  real :: vhtot(SZI_(G), SZJQ_(G))  ! The vertical sum of vhD, in m3 s-1.
  real, dimension(SZIQ_(G)) :: &
    T_u, S_u, &   ! Temperature, salinity, and pressure on the interface at
    pres_u        ! the u-point in the horizontal.
  real, dimension(SZI_(G)) :: &
    T_v, S_v, &   ! Temperature, salinity, and pressure on the interface at
    pres_v        ! the v-point in the horizontal.
  real :: Work_u(SZIQ_(G), SZJ_(G)) ! The work being done by the thickness
  real :: Work_v(SZI_(G), SZJQ_(G)) ! diffusion integrated over a cell, in W.
  real :: Work_h                    ! The work averaged over an h-cell in W m-2.
  real :: I4dt                      ! 1 / 4 dt
  real :: drdiA, drdiB  ! Along layer zonal- and meridional- potential density
  real :: drdjA, drdjB  ! gradients in the layers above (A) and below(B) the
                        ! interface times the grid spacing, in kg m-3.
  real :: drdkL, drdkR  ! Vertical density differences across an interface,
                        ! in kg m-3.
  real :: hg2A, hg2B, hg2L, hg2R
  real :: haA, haB, haL, haR
  real :: dzaL, dzaR
  real :: wtA, wtB, wtL, wtR
  real :: drdx, drdy, drdz  ! Zonal, meridional, and vertical density gradients,
                            ! in units of kg m-4.
  real :: Sfn_est       ! Two preliminary estimates (before limiting) of the
  real :: Sfn_unlim     ! overturning streamfunction, both in m3 s-1.
  real :: Sfn           ! The overturning streamfunction, in m3 s-1.
  real :: Sfn_safe      ! The streamfunction that goes linearly back to 0 at the
                        ! top.  This is a good thing to use when the slope is
                        ! so large as to be meaningless.
  real :: Slope         ! The slope of density surfaces, calculated in a way
                        ! that is always between -1 and 1.
  real :: mag_grad2     ! The squared magnitude of the 3-d density gradient, in kg2 m-8.
  real :: slope2_Ratio  ! The ratio of the slope squared to slope_max squared.
  real :: I_slope_max2  ! The inverse of slope_max squared, nondimensional.
  real :: h_neglect     ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected, in H.
  real :: h_neglect2    ! h_neglect^2, in H2.
  real :: dz_neglect    ! A thickness in m that is so small it is usually lost
                        ! in roundoff and can be neglected, in m.
  real :: G_scale       ! The gravitational accerlation times the conversion
                        ! factor from thickness to m, in m s-2 or m4 s-2 kg-1.
  logical :: use_EOS    ! If true, density is calculated from T & S using an
                        ! equation of state.
  logical :: find_work  ! If true, find the change in energy due to the fluxes.
  integer :: nk_linear  ! The number of layers over which the streamfunction
                        ! goes to 0.

! Diagnostics that should be eliminated altogether later...
 ! real, dimension(SZIQ_(G), SZJ_(G), SZK_(G)+1) :: sfn_x, sfn_slope_x
 ! real, dimension(SZI_(G), SZJQ_(G), SZK_(G)+1) :: sfn_y, sfn_slope_y

  integer :: is, ie, js, je, nz, Isdq
  integer :: i, j, k
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke ; Isdq = G%Isdq

  I4dt = 0.25 / dt
  I_slope_max2 = 1.0 / (CS%slope_max**2)
  G_scale = G%g_Earth * G%H_to_m
  h_neglect = G%H_subroundoff ; h_neglect2 = h_neglect**2
  dz_neglect = G%H_subroundoff*G%H_to_m

  use_EOS = associated(tv%eqn_of_state)

  nk_linear = max(G%nkml, 1)

  find_work = .false.
  if (ASSOCIATED(MEKE)) find_work = ASSOCIATED(MEKE%GM_src)
  find_work = (ASSOCIATED(CS%GMwork) .or. find_work)

  if (use_EOS) then
    call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth, dt, T, S, G, 1)
  endif

  ! Find the maximum and minimum permitted streamfunction.
  do j=js-1,je+1 ; do i=is-1,ie+1
    h_avail_rsum(i,j,1) = 0.0
    pres(i,j,1) = 0.0  ! ### This should be atmospheric pressure.

    h_avail(i,j,1) = max(I4dt*G%DXDYh(i,j)*(h(i,j,1)-G%Angstrom),0.0)
    h_avail_rsum(i,j,2) = h_avail(i,j,1)
    h_frac(i,j,1) = 1.0
    pres(i,j,2) = pres(i,j,1) + G%H_to_Pa*h(i,j,1)
  enddo ; enddo
  do k=2,nz
    do j=js-1,je+1 ; do i=is-1,ie+1
      h_avail(i,j,k) = max(I4dt*G%DXDYh(i,j)*(h(i,j,k)-G%Angstrom),0.0)
      h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
      h_frac(i,j,k) = 0.0 ; if (h_avail(i,j,k) > 0.0) &
        h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
      pres(i,j,K+1) = pres(i,j,K) + G%H_to_Pa*h(i,j,k)
    enddo ; enddo
  enddo

  do j=js,je ; do I=is-1,ie
    uhtot(I,j) = 0.0 ; Work_u(I,j) = 0.0
    ! sfn_x(I,j,1) = 0.0 ; sfn_slope_x(I,j,1) = 0.0
    ! sfn_x(I,j,nz+1) = 0.0 ; sfn_slope_x(I,j,nz+1) = 0.0
  enddo ; enddo
  do J=js-1,je ; do i=is,ie
    vhtot(i,J) = 0.0 ; Work_v(i,J) = 0.0
    ! sfn_y(i,J,1) = 0.0 ; sfn_slope_y(i,J,1) = 0.0
    ! sfn_y(i,J,nz+1) = 0.0 ; sfn_slope_y(i,J,nz+1) = 0.0
  enddo ; enddo

  do K=nz,2,-1
    if (find_work .and. .not.(use_EOS)) then
      drdiA = 0.0 ; drdiB = 0.0 ; drdjA = 0.0 ; drdjB = 0.0
!       drdkL = G%g_prime(k) ; drdkR = G%g_prime(k)
      drdkL = G%Rlay(k)-G%Rlay(k-1) ; drdkR = G%Rlay(k)-G%Rlay(k-1)
    endif

    ! Calculate the zonal fluxes and gradients.
    if (use_EOS .and. ((k > nk_linear) .or. find_work)) then
      do j=js,je
        do I=is-1,ie
          pres_u(I) = 0.5*(pres(i,j,K) + pres(i+1,j,K))
          T_u(I) = 0.25*((T(i,j,k) + T(i+1,j,k)) + (T(i,j,k-1) + T(i+1,j,k-1)))
          S_u(I) = 0.25*((S(i,j,k) + S(i+1,j,k)) + (S(i,j,k-1) + S(i+1,j,k-1)))
        enddo
        call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u(:,j), &
                     drho_dS_u(:,j), (is-Isdq+1)-1, ie-is+2, tv%eqn_of_state)
      enddo
    endif

    do j=js,je ; do I=is-1,ie
      if (use_EOS .and. ((k > nk_linear) .or. find_work)) then
        ! Estimate the horizontal density gradients along layers.
        drdiA = drho_dT_u(I,j) * (T(i+1,j,k-1)-T(i,j,k-1)) + &
                drho_dS_u(I,j) * (S(i+1,j,k-1)-S(i,j,k-1))
        drdiB = drho_dT_u(I,j) * (T(i+1,j,k)-T(i,j,k)) + &
                drho_dS_u(I,j) * (S(i+1,j,k)-S(i,j,k))

        ! Estimate the vertical density gradients times the grid spacing.
        drdkL = (drho_dT_u(I,j) * (T(i,j,k)-T(i,j,k-1)) + &
                 drho_dS_u(I,j) * (S(i,j,k)-S(i,j,k-1)))
        drdkR = (drho_dT_u(I,j) * (T(i+1,j,k)-T(i+1,j,k-1)) + &
                 drho_dS_u(I,j) * (S(i+1,j,k)-S(i+1,j,k-1)))
      endif


      if (k > nk_linear) then
        if (use_EOS) then
          hg2A = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
          hg2B = h(i,j,k)*h(i+1,j,k) + h_neglect2
          hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
          hg2R = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
          haA = 0.5*(h(i,j,k-1) + h(i+1,j,k-1))
          haB = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
          haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
          haR = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
          if (G%Boussinesq) then
            dzaL = haL * G%H_to_m ; dzaR = haR * G%H_to_m
          else
            dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect
            dzaR = 0.5*(e(i+1,j,K-1) - e(i+1,j,K+1)) + dz_neglect
          endif
          ! Use the harmonic mean thicknesses to weight the horizontal gradients.
          ! These unnormalized weights have been rearranged to minimize divisions.
          wtA = hg2A*haB ; wtB = hg2B*haA
          wtL = hg2L*(haR*dzaR) ; wtR = hg2R*(haL*dzaL)

          drdz = (wtL * drdkL + wtR * drdkR) / (dzaL*wtL + dzaR*wtR)
          ! The expression for drdz above is mathematically equivalent to:
          !   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
          !          ((hg2L/haL) + (hg2R/haR))
          ! This is the gradient of density along geopotentials.
          drdx = ((wtA * drdiA + wtB * drdiB) / (wtA + wtB) - &
                  drdz * (e(i,j,K)-e(i+1,j,K))) * G%IDXu(I,j)

          ! This estimate of slope is accurate for small slopes, but bounded
          ! to be between -1 and 1.
          mag_grad2 = drdx**2 + drdz**2
          if (mag_grad2 > 0.0) then
            Slope = drdx / sqrt(mag_grad2)
            slope2_Ratio = Slope**2 * I_slope_max2
          else ! Just in case mag_grad2 = 0 ever.
            Slope = 0.0
            slope2_Ratio = 1.0e20  ! Force the use of the safe streamfunction.
          endif

          if (present(int_slope_u)) then
            Slope = (1.0 - int_slope_u(I,j,K)) * Slope + &
                    int_slope_u(I,j,K) * ((e(i+1,j,K)-e(i,j,K)) * G%IDXu(I,j))
            slope2_Ratio = (1.0 - int_slope_u(I,j,K)) * slope2_Ratio
          endif

          ! Estimate the streamfunction at each interface.
          Sfn_unlim = -((KH_u(I,j,K)*G%dy_u(I,j))*Slope) * G%m_to_H
          if (uhtot(I,j) <= 0.0) then
            ! The transport that must balance the transport below is positive.
            Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i,j,k))
          else !  (uhtot(I,j) > 0.0)
            Sfn_safe = uhtot(I,j) * (1.0 - h_frac(i+1,j,k))
          endif

          ! Avoid moving dense water upslope from below the level of
          ! the bottom on the receiving side.
          if (Sfn_unlim > 0.0) then ! The flow below this interface is positive.
            if (e(i,j,K) < e(i+1,j,nz+1)) then
              Sfn_unlim = 0.0 ! This is not uhtot, because it may compensate for
                              ! deeper flow in very unusual cases.
            elseif (e(i+1,j,nz+1) > e(i,j,K+1)) then
              ! Scale the transport with the fraction of the donor layer above
              ! the bottom on the receiving side.
              Sfn_unlim = Sfn_unlim * ((e(i,j,K) - e(i+1,j,nz+1)) / &
                                       ((e(i,j,K) - e(i,j,K+1)) + dz_neglect))
            endif
          else
            if (e(i+1,j,K) < e(i,j,nz+1)) then ; Sfn_unlim = 0.0
            elseif (e(i,j,nz+1) > e(i+1,j,K+1)) then
              Sfn_unlim = Sfn_unlim * ((e(i+1,j,K) - e(i,j,nz+1)) / &
                                     ((e(i+1,j,K) - e(i+1,j,K+1)) + dz_neglect))
            endif
          endif

          Sfn_est = (Sfn_unlim + slope2_Ratio*Sfn_safe) / (1.0 + slope2_Ratio)
        else  ! With .not.use_EOS, the layers are constant density.
          Sfn_est = ((KH_u(I,j,K)*G%dy_u(I,j)) * &
                     ((e(i,j,K)-e(i+1,j,K))*G%IDXu(I,j))) * G%m_to_H
        endif

        ! Make sure that there is enough mass above to allow the streamfunction
        ! to satisfy the boundary condition of 0 at the surface.
        Sfn = min(max(Sfn_est, -h_avail_rsum(i,j,K)), h_avail_rsum(i+1,j,K))

        ! The actual transport is limited by the mass available in the two
        ! neighboring grid cells.
        uhD(I,j,k) = max(min((Sfn - uhtot(I,j)), h_avail(i,j,k)), &
                         -h_avail(i+1,j,k))

 !       sfn_x(I,j,K) = max(min(Sfn, uhtot(I,j)+h_avail(i,j,k)), &
 !                          uhtot(I,j)-h_avail(i+1,j,K))
 !       sfn_slope_x(I,j,K) = max(uhtot(I,j)-h_avail(i+1,j,k), &
 !                                min(uhtot(I,j)+h_avail(i,j,k), &
 !             min(h_avail_rsum(i+1,j,K), max(-h_avail_rsum(i,j,K), &
 !             (KH_u(I,j,K)*G%dy_u(I,j)) * ((e(i,j,K)-e(i+1,j,K))*G%IDXu(I,j)) )) ))
      else ! k <= nk_linear
        ! Balance the deeper flow with a return flow uniformly distributed
        ! though the remaining near-surface layers.  This is the same as
        ! using Sfn_safe above.  There is no need to apply the limiters in
        ! this case.
        if (uhtot(I,j) <= 0.0) then
          uhD(I,j,k) = -uhtot(I,j) * h_frac(i,j,k)
        else !  (uhtot(I,j) > 0.0)
          uhD(I,j,k) = -uhtot(I,j) * h_frac(i+1,j,k)
        endif

 !       sfn_x(I,j,K) = sfn_x(I,j,K+1) + uhD(I,j,k)
 !       if (sfn_slope_x(I,j,K+1) <= 0.0) then
 !         sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i,j,k))
 !       else
 !         sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i+1,j,k))
 !       endif
      endif

      uhtot(I,j) = uhtot(I,j) + uhD(I,j,k)

      if (find_work) then
        !   This is the energy tendency based on the original profiles, and does
        ! not include any nonlinear terms due to a finite time step (which would
        ! involve interactions between the fluxes through the different faces.
        !   A second order centered estimate is used for the density transfered
        ! between water columns.

        Work_u(I,j) = Work_u(I,j) + G_scale * &
          ( uhtot(I,j) * (drdkR * e(i+1,j,K) - drdkL * e(i,j,K)) - &
            (uhD(I,j,K) * drdiB) * 0.25 * &
            ((e(i,j,K) + e(i,j,K+1)) + (e(i+1,j,K) + e(i+1,j,K+1))) )
      endif

    enddo ; enddo

    ! Calculate the meridional fluxes and gradients.
    if (use_EOS .and. ((k > nk_linear) .or. find_work)) then
      do J=js-1,je
        do i=is,ie
          pres_v(i) = 0.5*(pres(i,j,K) + pres(i,j+1,K))
          T_v(i) = 0.25*((T(i,j,k) + T(i,j+1,k)) + (T(i,j,k-1) + T(i,j+1,k-1)))
          S_v(i) = 0.25*((S(i,j,k) + S(i,j+1,k)) + (S(i,j,k-1) + S(i,j+1,k-1)))
        enddo
        call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v(:,J), &
                     drho_dS_v(:,J), is, ie-is+1, tv%eqn_of_state)
      enddo
    endif
    do J=js-1,je ; do i=is,ie
      if (use_EOS .and. ((k > nk_linear) .or. find_work)) then
        ! Estimate the horizontal density gradients along layers.
        drdjA = drho_dT_v(i,J) * (T(i,j+1,k-1)-T(i,j,k-1)) + &
                drho_dS_v(i,J) * (S(i,j+1,k-1)-S(i,j,k-1))
        drdjB = drho_dT_v(i,J) * (T(i,j+1,k)-T(i,j,k)) + &
                drho_dS_v(i,J) * (S(i,j+1,k)-S(i,j,k))

        ! Estimate the vertical density gradients times the grid spacing.
        drdkL = (drho_dT_v(i,J) * (T(i,j,k)-T(i,j,k-1)) + &
                 drho_dS_v(i,J) * (S(i,j,k)-S(i,j,k-1)))
        drdkR = (drho_dT_v(i,J) * (T(i,j+1,k)-T(i,j+1,k-1)) + &
                 drho_dS_v(i,J) * (S(i,j+1,k)-S(i,j+1,k-1)))
      endif

      if (k > nk_linear) then
        if (use_EOS) then
          hg2A = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
          hg2B = h(i,j,k)*h(i,j+1,k) + h_neglect2
          hg2L = h(i,j,k-1)*h(i,j,k) + h_neglect2
          hg2R = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
          haA = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
          haB = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
          haL = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
          haR = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
          if (G%Boussinesq) then
            dzaL = haL * G%H_to_m ; dzaR = haR * G%H_to_m
          else
            dzaL = 0.5*(e(i,j,K-1) - e(i,j,K+1)) + dz_neglect
            dzaR = 0.5*(e(i,j+1,K-1) - e(i,j+1,K+1)) + dz_neglect
          endif
          ! Use the harmonic mean thicknesses to weight the horizontal gradients.
          ! These unnormalized weights have been rearranged to minimize divisions.
          wtA = hg2A*haB ; wtB = hg2B*haA
          wtL = hg2L*(haR*dzaR) ; wtR = hg2R*(haL*dzaL)

          drdz = (wtL * drdkL + wtR * drdkR) / (dzaL*wtL + dzaR*wtR)
          ! The expression for drdz above is mathematically equivalent to:
          !   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
          !          ((hg2L/haL) + (hg2R/haR))
          ! This is the gradient of density along geopotentials.
          drdy = ((wtA * drdjA + wtB * drdjB) / (wtA + wtB) - &
                  drdz * (e(i,j,K)-e(i,j+1,K))) * G%IDYv(i,J)

          ! This estimate of slope is accurate for small slopes, but bounded
          ! to be between -1 and 1.
          mag_grad2 = drdy**2 + drdz**2
          if (mag_grad2 > 0.0) then
            Slope = drdy / sqrt(mag_grad2)
            slope2_Ratio = Slope**2 * I_slope_max2
          else ! Just in case mag_grad2 = 0 ever.
            Slope = 0.0
            slope2_Ratio = 1.0e20  ! Force the use of the safe streamfunction.
          endif

          if (present(int_slope_v)) then
            Slope = (1.0 - int_slope_v(i,J,K)) * Slope + &
                    int_slope_v(i,J,K) * ((e(i,j+1,K)-e(i,j,K)) * G%IDYv(i,J))
            slope2_Ratio = (1.0 - int_slope_v(i,J,K)) * slope2_Ratio
          endif

          ! Estimate the streamfunction at each interface.
          Sfn_unlim = -((KH_v(i,J,K)*G%dx_v(i,J))*Slope) * G%m_to_H
          if (vhtot(i,J) <= 0.0) then
            Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j,k))
          else !  (vhtot(I,j) > 0.0)
            Sfn_safe = vhtot(i,J) * (1.0 - h_frac(i,j+1,k))
          endif

          ! Avoid moving dense water upslope from below the level of
          ! the bottom on the receiving side.
          if (Sfn_unlim > 0.0) then ! The flow below this interface is positive.
            if (e(i,j,K) < e(i,j+1,nz+1)) then
              Sfn_unlim = 0.0 ! This is not vhtot, because it may compensate for
                              ! deeper flow in very unusual cases.
            elseif (e(i,j+1,nz+1) > e(i,j,K+1)) then
              ! Scale the transport with the fraction of the donor layer above
              ! the bottom on the receiving side.
              Sfn_unlim = Sfn_unlim * ((e(i,j,K) - e(i,j+1,nz+1)) / &
                                       ((e(i,j,K) - e(i,j,K+1)) + dz_neglect))
            endif
          else
            if (e(i,j+1,K) < e(i,j,nz+1)) then ; Sfn_unlim = 0.0
            elseif (e(i,j,nz+1) > e(i,j+1,K+1)) then
              Sfn_unlim = Sfn_unlim * ((e(i,j+1,K) - e(i,j,nz+1)) / &
                                     ((e(i,j+1,K) - e(i,j+1,K+1)) + dz_neglect))
            endif
          endif

          ! Estimate the streamfunction at each interface.
          Sfn_est = (Sfn_unlim + slope2_Ratio*Sfn_safe) / (1.0 + slope2_Ratio)
        else      ! With .not.use_EOS, the layers are constant density.
          Sfn_est = ((KH_v(i,J,K)*G%dx_v(i,J)) * &
                     ((e(i,j,K)-e(i,j+1,K))*G%IDYv(i,J))) * G%m_to_H
        endif

        ! Make sure that there is enough mass above to allow the streamfunction
        ! to satisfy the boundary condition of 0 at the surface.
        Sfn = min(max(Sfn_est, -h_avail_rsum(i,j,K)), h_avail_rsum(i,j+1,K))

        ! The actual transport is limited by the mass available in the two
        ! neighboring grid cells.
        vhD(i,J,k) = max(min((Sfn - vhtot(i,J)), h_avail(i,j,k)), &
                         -h_avail(i,j+1,k))

  !      sfn_y(i,J,K) = max(min(Sfn, vhtot(i,J)+h_avail(i,j,k)), &
  !                         vhtot(i,J)-h_avail(i,j+1,k))
  !      sfn_slope_y(i,J,K) = max(vhtot(i,J)-h_avail(i,j+1,k), &
  !                               min(vhtot(i,J)+h_avail(i,j,k), &
  !            min(h_avail_rsum(i,j+1,K), max(-h_avail_rsum(i,j,K), &
  !            (KH_v(i,J,K)*G%dx_v(i,J)) * ((e(i,j,K)-e(i,j+1,K))*G%IDYv(i,J)) )) ))
      else  ! k <= nk_linear
        ! Balance the deeper flow with a return flow uniformly distributed
        ! though the remaining near-surface layers.
        if (vhtot(i,J) <= 0.0) then
          vhD(i,J,k) = -vhtot(i,J) * h_frac(i,j,k)
        else !  (vhtot(i,J) > 0.0)
          vhD(i,J,k) = -vhtot(i,J) * h_frac(i,j+1,k)
        endif

  !      sfn_y(i,J,K) = sfn_y(i,J,K+1) + vhD(i,J,k)
  !      if (sfn_slope_y(i,J,K+1) <= 0.0) then
  !        sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j,k))
  !      else
  !        sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j+1,k))
  !      endif
      endif

      vhtot(i,J) = vhtot(i,J)  + vhD(i,J,k)

      if (find_work) then
        !   This is the energy tendency based on the original profiles, and does
        ! not include any nonlinear terms due to a finite time step (which would
        ! involve interactions between the fluxes through the different faces.
        !   A second order centered estimate is used for the density transfered
        ! between water columns.

        Work_v(i,J) = Work_v(i,J) + G_scale * &
          ( vhtot(i,J) * (drdkR * e(i,j+1,K) - drdkL * e(i,j,K)) - &
           (vhD(i,J,K) * drdjB) * 0.25 * &
           ((e(i,j,K) + e(i,j,K+1)) + (e(i,j+1,K) + e(i+1,j,K+1))) )
      endif
    enddo ; enddo
  enddo ! k-loop

  ! In layer 1, enforce the boundary conditions that Sfn(z=0) = 0.0
  if (.not.find_work .or. .not.(use_EOS)) then
    do j=js,je ; do I=is-1,ie ; uhD(I,j,1) = -uhtot(I,j) ; enddo ; enddo
    do J=js-1,je ; do i=is,ie ; vhD(i,J,1) = -vhtot(i,J) ; enddo ; enddo
  else
    if (use_EOS) then ; do j=js,je
      do I=is-1,ie
        pres_u(I) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
        T_u(I) = 0.5*(T(i,j,1) + T(i+1,j,1))
        S_u(I) = 0.5*(S(i,j,1) + S(i+1,j,1))
      enddo
      call calculate_density_derivs(T_u, S_u, pres_u, drho_dT_u(:,j), &
                   drho_dS_u(:,j), (is-Isdq+1)-1, ie-is+2, tv%eqn_of_state)
    enddo ; endif
    do j=js,je ; do I=is-1,ie
      uhD(I,j,1) = -uhtot(I,j)

      if (use_EOS) then
        drdiB = drho_dT_u(I,j) * (T(i+1,j,1)-T(i,j,1)) + &
                drho_dS_u(I,j) * (S(i+1,j,1)-S(i,j,1))
      endif
      Work_u(I,j) = Work_u(I,j) + G_scale * ( (uhD(I,j,1) * drdiB) * 0.25 * &
          ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )

    enddo ; enddo

    if (use_EOS) then ; do J=js-1,je
      do i=is,ie
        pres_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
        T_v(i) = 0.5*(T(i,j,1) + T(i,j+1,1))
        S_v(i) = 0.5*(S(i,j,1) + S(i,j+1,1))
      enddo
      call calculate_density_derivs(T_v, S_v, pres_v, drho_dT_v(:,J), &
                   drho_dS_v(:,J), is, ie-is+1, tv%eqn_of_state)
    enddo ; endif
    do J=js-1,je ; do i=is,ie
      vhD(i,J,1) = -vhtot(i,J)

      if (use_EOS) then
        drdjB = drho_dT_v(i,J) * (T(i,j+1,1)-T(i,j,1)) + &
                drho_dS_v(i,J) * (S(i,j+1,1)-S(i,j,1))
      endif
      Work_v(i,J) = Work_v(i,J) - G_scale * ( (vhD(i,J,1) * drdjB) * 0.25 * &
          ((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
    enddo ; enddo
  endif

  if (find_work) then ; do j=js,je ; do i=is,ie
    ! Note that the units of Work_v and Work_u are W, while Work_h is W m-2.
    Work_h = 0.5 * G%IDXDYh(i,j) * &
      ((Work_u(I-1,j) + Work_u(I,j)) + (Work_v(i,J-1) + Work_v(i,J)))
    if (ASSOCIATED(CS%GMwork)) CS%GMwork(i,j) = Work_h
    if (ASSOCIATED(MEKE)) then ; if (ASSOCIATED(MEKE%GM_src)) then
      MEKE%GM_src(i,j) = MEKE%GM_src(i,j) + Work_h
    endif ; endif
  enddo ; enddo ; endif

 ! if (CS%id_sfn_x > 0) call post_data(CS%id_sfn_x, sfn_x, CS%diag)
 ! if (CS%id_sfn_y > 0) call post_data(CS%id_sfn_y, sfn_y, CS%diag)
 ! if (CS%id_sfn_slope_x > 0) call post_data(CS%id_sfn_slope_x, sfn_slope_x, CS%diag)
 ! if (CS%id_sfn_slope_y > 0) call post_data(CS%id_sfn_slope_y, sfn_slope_y, CS%diag)

end subroutine thickness_diffuse_full

subroutine add_detangling_Kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, CS, &
                             int_slope_u, int_slope_v)
  real, dimension(NXMEM_,NYMEM_,NZ_),    intent(in)    :: h
  real, dimension(NXMEM_,NYMEM_,NZp1_),  intent(in)    :: e
  real, dimension(NXMEMQ_,NYMEM_,NZp1_), intent(inout) :: Kh_u
  real, dimension(NXMEM_,NYMEMQ_,NZp1_), intent(inout) :: Kh_v
  real, dimension(NXMEMQ_,NYMEM_),       intent(in)    :: Kh_u_CFL
  real, dimension(NXMEM_,NYMEMQ_),       intent(in)    :: Kh_v_CFL
  type(thermo_var_ptrs),                 intent(in)    :: tv
  real,                                  intent(in)    :: dt
  type(ocean_grid_type),                 intent(in)    :: G
  type(thickness_diffuse_CS),            pointer       :: CS
  real, dimension(NXMEMQ_,NYMEM_,NZp1_), intent(inout) :: int_slope_u
  real, dimension(NXMEM_,NYMEMQ_,NZp1_), intent(inout) :: int_slope_v
! Arguments: h - Layer thickness, in m.
!  (in)      e - Interface heights relative to mean sea level, in m.
!  (inout)   Kh_u - Thickness diffusivity on interfaces at u points, in m2 s-1.
!  (inout)   Kh_v - Thickness diffusivity on interfaces at v points, in m2 s-1.
!  (in)      Kh_u_CFL - Maximum stable thickness diffusivity at u points, in m2 s-1.
!  (in)      Kh_v_CFL - Maximum stable thickness diffusivity at v points, in m2 s-1.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      dt - Time increment in s.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 thickness_diffuse_init.
!  (inout)   int_slope_u - Ratios that determine how much of the isopycnal
!                 slopes are taken directly from the interface slopes without
!                 consideration of density gradients.
!  (inout)   int_slope_v - Ratios that determine how much of the isopycnal
!                 slopes are taken directly from the interface slopes without
!                 consideration of density gradients.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    de_top     ! The distances between the top of a layer and the top of the
               ! region where the detangling is applied, in H.
  real, dimension(SZIQ_(G),SZJ_(G),SZK_(G)) :: &
    Kh_lay_u   ! The tentative interface height diffusivity for each layer at
               ! u points, in m2 s-1.
  real, dimension(SZI_(G),SZJQ_(G),SZK_(G)) :: &
    Kh_lay_v   ! The tentative interface height diffusivity for each layer at
               ! v points, in m2 s-1.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    de_bot     ! The distances from the bottom of the region where the
               ! detangling is applied, in H.
  real :: h1, h2    ! The thinner and thicker surrounding thicknesses, in H,
                    ! with the thinner modified near the boundaries to mask out
                    ! thickness variations due to topography, etc.
  real :: jag_Rat   ! The nondimensional jaggedness ratio for a layer, going
                    ! from 0 (smooth) to 1 (jagged).  This is the difference
                    ! between the arithmetic and harmonic mean thicknesses
                    ! normalized by the arithmetic mean thickness.
  real :: Kh_scale  ! A ratio by which Kh_u_CFL is scaled for maximally jagged
                    ! layers, nondim.
  real :: Kh_det    ! The detangling diffusivity, in m2 s-1.
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected, in H.

  real :: I_sl      ! The absolute value of the larger in magnitude of the slopes
                    ! above and below.
  real :: Rsl       ! The ratio of the smaller magnitude slope to the larger
                    ! magnitude one, ND. 0 <= Rsl <1.
  real :: IRsl      ! The (limited) inverse of Rsl, ND. 1 < IRsl <= 1e9.
  real :: dH        ! The thickness gradient divided by the damping timescale
                    ! and the ratio of the face length to the adjacent cell
                    ! areas for comparability with the diffusivities, in m2 s-1.
  real :: adH       ! The absolute value of dH, in m2 s-1.
  real :: sign      ! 1 or -1, with the same sign as the layer thickness gradient. 
  real :: sl_K      ! The sign-corrected slope of the interface above, ND.
  real :: sl_Kp1    ! The sign-corrected slope of the interface below, ND.
  real :: I_sl_K    ! The (limited) inverse of sl_K, ND.
  real :: I_sl_Kp1  ! The (limited) inverse of sl_Kp1, ND.
  real :: I_4t      ! A quarter of inverse of the damping timescale, in s-1.
  real :: Fn_R      ! A function of Rsl, such that Rsl < Fn_R < 1.
  real :: denom, I_denom ! A denominator and its inverse, various units.
  real :: Kh_min    ! A local floor on the diffusivity, in m2 s-1.
  real :: Kh_max    ! A local ceiling on the diffusivity, in m2 s-1.
  real :: wt1, wt2  ! Nondimensional weights.
  !   Variables used only in testing code.
  ! real, dimension(SZK_(G)) :: uh_here
  ! real, dimension(SZK_(G)+1) :: Sfn
  real :: dKh       ! An increment in the diffusivity, in m2 s-1.

  real, dimension(SZIQ_(G),SZK_(G)+1) :: &
    Kh_bg, &        ! The background (floor) value of Kh, in m2 s-1.
    Kh, &           ! The tentative value of Kh, in m2 s-1.
    Kh_detangle, &  ! The detangling diffusivity that could be used, in m2 s-1.
    Kh_min_max_p, & ! The smallest ceiling that can be placed on Kh(I,K)
                    ! based on the value of Kh(I,K+1), in m2 s-1.
    Kh_min_max_m, & ! The smallest ceiling that can be placed on Kh(I,K)
                    ! based on the value of Kh(I,K-1), in m2 s-1.
    ! The following are variables that define the relationships between
    ! successive values of Kh.
    ! Search for Kh that satisfy...
    !    Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
    !    Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
    !    Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
    !    Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)
    Kh_min_m , &   ! See above, ND.
    Kh0_min_m , &  ! See above, in m2 s-1.
    Kh_max_m , &   ! See above, ND.
    Kh0_max_m, &   ! See above, in m2 s-1.
    Kh_min_p , &   ! See above, ND.
    Kh0_min_p , &  ! See above, in m2 s-1.
    Kh_max_p , &   ! See above, ND.
    Kh0_max_p      ! See above, in m2 s-1.
  real, dimension(SZIQ_(G)) :: &
    Kh_max_max  ! The maximum diffusivity permitted in a column.
  logical, dimension(SZIQ_(G)) :: &
    do_i        ! If true, work on a column.
  integer :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k_top
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke

  k_top = G%nk_rho_varies + 1
  h_neglect = G%H_subroundoff
  !   The 0.5 is because we are not using uniform weightings, but are
  ! distributing the diffusivities more effectively (with wt1 & wt2), but this
  ! means that the additions to a single interface can be up to twice as large.
  Kh_scale = 0.5
  if (CS%detangle_time > dt) Kh_scale = 0.5 * dt / CS%detangle_time

  do j=js-1,je+1 ; do i=is-1,ie+1
    de_top(i,j,k_top) = 0.0 ; de_bot(i,j) = 0.0
  enddo ; enddo
  do k=k_top+1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
    de_top(i,j,k) = de_top(i,j,k-1) + h(i,j,k-1)
  enddo ; enddo ; enddo

  do j=js,je ; do I=is-1,ie
    Kh_lay_u(I,j,nz) = 0.0 ; Kh_lay_u(I,j,k_top) = 0.0
  enddo ; enddo
  do J=js-1,je ; do i=is,ie
    Kh_lay_v(i,J,nz) = 0.0 ; Kh_lay_v(i,J,k_top) = 0.0
  enddo ; enddo

  do k=nz-1,k_top+1,-1
    ! Find the diffusivities associated with each layer.
    do j=js-1,je+1 ; do i=is-1,ie+1
      de_bot(i,j) = de_bot(i,j) + h(i,j,k+1)
    enddo ; enddo

    do j=js,je ; do I=is-1,ie ; if (G%umask(I,j) > 0.0) then
      if (h(i,j,k) > h(i+1,j,k)) then
        h2 = h(i,j,k)
        h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) )
      else
        h2 = h(i+1,j,k)
        h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
      endif
      jag_Rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
      Kh_lay_u(I,j,k) = (Kh_scale * Kh_u_CFL(I,j)) * jag_Rat**2
    endif ; enddo ; enddo

    do J=js-1,je ; do i=is,ie ; if (G%vmask(i,J) > 0.0) then
      if (h(i,j,k) > h(i,j+1,k)) then
        h2 = h(i,j,k)
        h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) )
      else
        h2 = h(i,j+1,k)
        h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
      endif
      jag_Rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
      Kh_lay_v(i,J,k) = (Kh_scale * Kh_v_CFL(i,J)) * jag_Rat**2
    endif ; enddo ; enddo
  enddo

  ! Limit the diffusivities
  
  I_4t = Kh_scale / (4.0*dt)

  do n=1,2
    if (n==1) then ; jsh = js ; ish = is-1
    else ; jsh = js-1 ; ish = is ; endif
  
    do j=jsh,je

      ! First, populate the diffusivities
      if (n==1) then ! This is a u-column.
        do i=ish,ie
          do_i(I) = (G%umask(I,j) > 0.0)
          Kh_max_max(I) = Kh_u_CFL(I,j)
        enddo
        do K=1,nz+1 ; do i=ish,ie
          Kh_bg(I,K) = Kh_u(I,j,K) ; Kh(I,K) = Kh_bg(I,K)
          Kh_min_max_p(I,K) = Kh_bg(I,K) ; Kh_min_max_m(I,K) = Kh_bg(I,K)
          Kh_detangle(I,K) = 0.0
        enddo ; enddo
      else ! This is a v-column.
        do i=ish,ie
          do_i(i) = (G%vmask(i,J) > 0.0) ; Kh_max_max(I) = Kh_v_CFL(i,J)
        enddo
        do K=1,nz+1 ; do i=ish,ie
          Kh_bg(I,K) = Kh_v(I,j,K) ; Kh(I,K) = Kh_bg(I,K)
          Kh_min_max_p(I,K) = Kh_bg(I,K) ; Kh_min_max_m(I,K) = Kh_bg(I,K)
          Kh_detangle(I,K) = 0.0
        enddo ; enddo
      endif

      ! Determine the limits on the diffusivities.
      do k=k_top,nz ; do i=ish,ie ; if (do_i(i)) then
        if (n==1) then ! This is a u-column.
          dH = 0.0
          denom = ((G%IDXDYh(i+1,j) + G%IDXDYh(i,j))*G%dy_u(I,j))
          !   This expression uses differences in e in place of h for better
          ! consistency with the slopes.
          if (denom > 0.0) &
            dH = I_4t * ((e(i+1,j,K) - e(i+1,j,K+1)) - &
                         (e(i,j,K) - e(i,j,K+1))) / denom
           ! dH = I_4t * (h(i+1,j,k) - h(i,j,k)) / denom

          adH = abs(dH)
          sign = 1.0 ; if (dH < 0) sign = -1.0
          sl_K = sign * (e(i+1,j,K)-e(i,j,K)) * G%IDXu(I,j)
          sl_Kp1 = sign * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IDXu(I,j)

          ! Add the incremental diffusivites to the surrounding interfaces.
          ! Adding more to the more steeply sloping layers (as below) makes
          ! the diffusivities more than twice as effective.
          denom = (sl_K**2 + sl_Kp1**2)
          wt1 = 0.5 ; wt2 = 0.5
          if (denom > 0.0) then
            wt1 = sl_K**2 / denom ; wt2 = sl_Kp1**2 / denom
          endif
          Kh_detangle(I,K) = Kh_detangle(I,K) + wt1*Kh_lay_u(I,j,k)
          Kh_detangle(I,K+1) = Kh_detangle(I,K+1) + wt2*Kh_lay_u(I,j,k)
        else ! This is a v-column.
          dH = 0.0
          denom = ((G%IDXDYh(i,j+1) + G%IDXDYh(i,j))*G%dx_v(I,j))
          if (denom > 0.0) &
            dH = I_4t * ((e(i,j+1,K) - e(i,j+1,K+1)) - &
                         (e(i,j,K) - e(i,j,K+1))) / denom
           ! dH = I_4t * (h(i,j+1,k) - h(i,j,k)) / denom

          adH = abs(dH)
          sign = 1.0 ; if (dH < 0) sign = -1.0
          sl_K = sign * (e(i,j+1,K)-e(i,j,K)) * G%IDYv(i,J)
          sl_Kp1 = sign * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IDYv(i,J)

          ! Add the incremental diffusviites to the surrounding interfaces.
          ! Adding more to the more steeply sloping layers (as below) makes
          ! the diffusivities more than twice as effective.
          denom = (sl_K**2 + sl_Kp1**2)
          wt1 = 0.5 ; wt2 = 0.5
          if (denom > 0.0) then
            wt1 = sl_K**2 / denom ; wt2 = sl_Kp1**2 / denom
          endif
          Kh_detangle(I,K) = Kh_detangle(I,K) + wt1*Kh_lay_v(i,J,k)
          Kh_detangle(I,K+1) = Kh_detangle(I,K+1) + wt2*Kh_lay_v(i,J,k)
        endif

        if (adH == 0.0) then
          Kh_min_m(I,K+1) = 1.0 ; Kh0_min_m(I,K+1) = 0.0
          Kh_max_m(I,K+1) = 1.0 ; Kh0_max_m(I,K+1) = 0.0
          Kh_min_p(I,K) = 1.0 ; Kh0_min_p(I,K) = 0.0
          Kh_max_p(I,K) = 1.0 ; Kh0_max_p(I,K) = 0.0
        elseif (adH > 0.0) then
          if (sl_K <= sl_Kp1) then
            ! This case should only arise from nonlinearities in the equation of state.
            ! Treat it as though dedx(K) = dedx(K+1) & dH = 0.
            Kh_min_m(I,K+1) = 1.0 ; Kh0_min_m(I,K+1) = 0.0
            Kh_max_m(I,K+1) = 1.0 ; Kh0_max_m(I,K+1) = 0.0
            Kh_min_p(I,K) = 1.0 ; Kh0_min_p(I,K) = 0.0
            Kh_max_p(I,K) = 1.0 ; Kh0_max_p(I,K) = 0.0
          elseif (sl_K <= 0.0) then   ! Both slopes are opposite to dH
            I_sl = -1.0 / sl_Kp1
            Rsl = -sl_K * I_sl                            ! 0 <= Rsl < 1 
            IRsl = 1e9 ; if (Rsl > 1e-9) IRsl = 1.0/Rsl   ! 1 < IRsl <= 1e9

            Fn_R = Rsl
            if (Kh_max_max(I) > 0) &
              Fn_R = min(sqrt(Rsl), Rsl + (adH * I_sl) / Kh_max_max(I))

            Kh_min_m(I,K+1) = Fn_R ; Kh0_min_m(I,K+1) = 0.0
            Kh_max_m(I,K+1) = Rsl ; Kh0_max_m(I,K+1) = adH * I_sl
            Kh_min_p(I,K) = IRsl ; Kh0_min_p(I,K) = -adH * (I_sl*IRsl)
            Kh_max_p(I,K) = 1.0/(Fn_R + 1.0e-30) ; Kh0_max_p(I,K) = 0.0
          elseif (sl_Kp1 < 0.0) then  ! Opposite (nonzero) signs of slopes.
            I_sl_K = 1e18 ; if (sl_K > 1e-18) I_sl_K = 1.0 / sl_K
            I_sl_Kp1 = 1e18 ; if (-sl_Kp1 > 1e-18) I_sl_Kp1 = -1.0 / sl_Kp1

            Kh_min_m(I,K+1) = 0.0 ; Kh0_min_m(I,K+1) = 0.0
            Kh_max_m(I,K+1) = - sl_K*I_sl_Kp1 ; Kh0_max_m(I,K+1) = adH*I_sl_Kp1
            Kh_min_p(I,K) = 0.0 ; Kh0_min_p(I,K) = 0.0
            Kh_max_p(I,K) = sl_Kp1*I_sl_K ; Kh0_max_p(I,K) = adH*I_sl_K

            ! This limit does not use the slope weighting so that potentially
            ! sharp gradients in diffusivities are not forced to occur.
            Kh_max = adH / (sl_K - sl_Kp1)
            Kh_min_max_p(I,K) = max(Kh_min_max_p(I,K), Kh_max)
            Kh_min_max_m(I,K+1) = max(Kh_min_max_m(I,K+1), Kh_max)
          else ! Both slopes are of the same sign as dH.
            I_sl = 1.0 / sl_K
            Rsl = sl_Kp1 * I_sl                           ! 0 <= Rsl < 1 
            IRsl = 1e9 ; if (Rsl > 1e-9) IRsl = 1.0/Rsl   ! 1 < IRsl <= 1e9

            ! Rsl <= Fn_R <= 1
            Fn_R = Rsl
            if (Kh_max_max(I) > 0) &
              Fn_R = min(sqrt(Rsl), Rsl + (adH * I_sl) / Kh_max_max(I))

            Kh_min_m(I,K+1) = IRsl ; Kh0_min_m(I,K+1) = -adH * (I_sl*IRsl)
            Kh_max_m(I,K+1) = 1.0/(Fn_R + 1.0e-30) ; Kh0_max_m(I,K+1) = 0.0
            Kh_min_p(I,K) = Fn_R ; Kh0_min_p(I,K) = 0.0
            Kh_max_p(I,K) = Rsl ; Kh0_max_p(I,K) = adH * I_sl
          endif
        endif
      endif ; enddo ; enddo ! I-loop & k-loop

      do k=k_top,nz+1,nz+1-k_top ; do i=ish,ie ; if (do_i(i)) then
        ! The diffusivities at k_top and nz+1 are both fixed. 
        Kh_min_m(I,k) = 0.0 ; Kh0_min_m(I,k) = 0.0
        Kh_max_m(I,k) = 0.0 ; Kh0_max_m(I,k) = 0.0
        Kh_min_p(I,k) = 0.0 ; Kh0_min_p(I,k) = 0.0
        Kh_max_p(I,k) = 0.0 ; Kh0_max_p(I,k) = 0.0
        Kh_min_max_p(I,K) = Kh_bg(I,K)
        Kh_min_max_m(I,K) = Kh_bg(I,K)
      endif ; enddo ; enddo ! I-loop and k_top/nz+1 loop
 
      ! Search for Kh that satisfy...
      !    Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
      !    Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
      !    Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
      !    Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)

      ! Increase the diffusivies to satisfy the min constraints.
      ! All non-zero min constraints on one diffusivity are max constraints on another.
      do K=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
        Kh(I,K) = max(Kh_bg(I,K), Kh_detangle(I,K), &
                      min(Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K), Kh(I,K-1)))

        if (Kh0_max_m(I,K) > Kh_bg(I,K)) Kh(I,K) = min(Kh(I,K), Kh0_max_m(I,K))
        if (Kh0_max_p(I,K) > Kh_bg(I,K)) Kh(I,K) = min(Kh(I,K), Kh0_max_p(I,K))
      endif ; enddo ; enddo ! I-loop & k-loop
      ! This is still true... do i=ish,ie ; Kh(I,nz+1) = Kh_bg(I,nz+1) ; enddo
      do K=nz,k_top+1,-1 ; do i=ish,ie ; if (do_i(i)) then
        Kh(I,k) = max(Kh(I,K), min(Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K), Kh(I,K+1)))

        Kh_max = max(Kh_min_max_p(I,K), Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K))
        Kh(I,k) = min(Kh(I,k), Kh_max)
      endif ; enddo ; enddo ! I-loop & k-loop
      !  All non-zero min constraints on one diffusivity are max constraints on
      ! another layer, so the min constraints can now be discounted.

      ! Decrease the diffusivities to satisfy the max constraints.
        do K=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
          Kh_max = max(Kh_min_max_m(I,K), Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K))
          if (Kh(I,k) > Kh_max) Kh(I,k) = Kh_Max
        endif ; enddo ; enddo  ! i- and K-loops

      ! This code tests the solutions...
!     do i=ish,ie
!       Sfn(:) = 0.0 ; uh_here(:) = 0.0
!       do K=k_top,nz
!         if ((Kh(i,K) > Kh_bg(i,K)) .or. (Kh(i,K+1) > Kh_bg(i,K+1))) then
!           if (n==1) then ! u-point.
!             if ((h(i+1,j,k) - h(i,j,k)) * &
!                 ((e(i+1,j,K)-e(i+1,j,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
!               Sfn(K) = -Kh(i,K) * (e(i+1,j,K)-e(i,j,K)) * G%IDXu(I,j)
!               Sfn(K+1) = -Kh(i,K+1) * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IDXu(I,j)
!               uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dy_u(I,j)
!               if (abs(uh_here(k))*min(G%IDXDYh(i,j), G%IDXDYh(i+1,j)) > &
!                   (1e-10*G%m_to_H)) then
!                 if (uh_here(k) * (h(i+1,j,k) - h(i,j,k)) > 0.0) then
!                   call GOLD_error(WARNING, &
!                          "Corrective u-transport is up the thickness gradient.", .true.)
!                 endif
!                 if (((h(i,j,k) - 4.0*dt*G%IDXDYh(i,j)*uh_here(k)) - &
!                      (h(i+1,j,k) + 4.0*dt*G%IDXDYh(i+1,j)*uh_here(k))) * &
!                     (h(i,j,k) - h(i+1,j,k)) < 0.0) then
!                   call GOLD_error(WARNING, &
!                          "Corrective u-transport is too large.", .true.)
!                 endif
!               endif
!             endif
!           else ! v-point
!             if ((h(i,j+1,k) - h(i,j,k)) * &
!                 ((e(i,j+1,K)-e(i,j+1,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
!               Sfn(K) = -Kh(i,K) * (e(i,j+1,K)-e(i,j,K)) * G%IDYv(i,J)
!               Sfn(K+1) = -Kh(i,K+1) * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IDYv(i,J)
!               uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dx_v(i,J)
!               if (abs(uh_here(K))*min(G%IDXDYh(i,j), G%IDXDYh(i,j+1)) > &
!                   (1e-10*G%m_to_H)) then
!                 if (uh_here(K) * (h(i,j+1,k) - h(i,j,k)) > 0.0) then
!                   call GOLD_error(WARNING, &
!                          "Corrective v-transport is up the thickness gradient.", .true.)
!                 endif
!                 if (((h(i,j,k) - 4.0*dt*G%IDXDYh(i,j)*uh_here(K)) - &
!                      (h(i,j+1,k) + 4.0*dt*G%IDXDYh(i,j+1)*uh_here(K))) * &
!                     (h(i,j,k) - h(i,j+1,k)) < 0.0) then
!                   call GOLD_error(WARNING, &
!                          "Corrective v-transport is too large.", .true.)
!                 endif
!               endif
!             endif
!           endif ! u- or v- selection.
!          !  de_dx(I,K) = (e(i+1,j,K)-e(i,j,K)) * G%IDXu(I,j)
!         endif
!       enddo
!     enddo

      if (n==1) then ! This is a u-column.
        do K=k_top+1,nz ; do i=ish,ie
          if (Kh(I,K) > Kh_u(I,j,K)) then
            dKh = (Kh(I,K) - Kh_u(I,j,K))
            int_slope_u(I,j,K) = dKh / Kh(I,K)
            Kh_u(I,j,K) = Kh(I,K)
          endif
        enddo ; enddo
      else ! This is a v-column.
        do K=k_top+1,nz ; do i=ish,ie
          if (Kh(i,K) > Kh_v(i,J,K)) then
            dKh = Kh(i,K) - Kh_v(i,J,K)
            int_slope_v(i,J,K) = dKh / Kh(i,K)
            Kh_v(i,J,K) = Kh(i,K)
          endif
        enddo ; enddo
      endif

    enddo ! j-loop
  enddo  ! n-loop over u- and v- directions.

end subroutine add_detangling_Kh

subroutine vert_fill_TS(h, T_in, S_in, kappa, dt, T_f, S_f, G, halo_here)
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(in)    :: h
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(in)    :: T_in
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(in)    :: S_in
  real,                                intent(in)    :: kappa
  real,                                intent(in)    :: dt
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(out)   :: T_f
  real, dimension(NXMEM_,NYMEM_,NZ_),  intent(out)   :: S_f
  type(ocean_grid_type),               intent(in)    :: G
  integer,                   optional, intent(in)    :: halo_here
!    This subroutine fills massless layers with sensible values of two
!*  tracer arrays (nominally temperature and salinity) by diffusing
!*  vertically with a (small?) constant diffusivity.

! Arguments: h - Layer thickness, in m or kg m-2.
!  (in)      T_in - The input temperature, in K.
!  (in)      S_in - The input salinity, in psu.
!  (in)      kappa - The diapycnal diffusivity, in m2 s-1.
!  (in)      dt - Time increment in s.
!  (out)     T_f - The filled temperature, in K.
!  (out)     S_f - The filled salinity, in psu.
!  (in)      G - The ocean's grid structure.
!  (in,opt)  halo_here - the number of halo points to work on, 0 by default.

  real :: ent(SZI_(G),SZK_(G)+1)   ! The diffusive entrainment (kappa*dt)/dz
                                   ! between layers in a timestep in m or kg m-2.
  real :: b1(SZI_(G)), d1(SZI_(G)) ! b1, c1, and d1 are variables used by the
  real :: c1(SZI_(G),SZK_(G))      ! tridiagonal solver.
  real :: kap_dt_x2                ! The product of 2*kappa*dt, converted to
                                   ! the same units as h, in m2 or kg2 m-4.
  real :: h0                       ! A negligible thickness, in m or kg m-2, to
                                   ! allow for zero thicknesses.
  integer :: i, j, k, is, ie, js, je, nz, halo

  halo=0 ; if (present(halo_here)) halo = max(halo_here,0)

  is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo
  nz = G%ke

  kap_dt_x2 = (2.0*kappa*dt)*G%m_to_H**2
  h0 = 1.0e-16*sqrt(kappa*dt)*G%m_to_H

  if (kap_dt_x2 <= 0.0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      T_f(i,j,k) = T_in(i,j,k) ; S_f(i,j,k) = S_in(i,j,k)
    enddo ; enddo ; enddo
  else
    do j=js,je
      do i=is,ie
        ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
        b1(i) = 1.0 / (h(i,j,1)+ent(i,2))
        d1(i) = b1(i) * h(i,j,1)
        T_f(i,j,1) = (b1(i)*h(i,j,1))*T_in(i,j,1)
        S_f(i,j,1) = (b1(i)*h(i,j,1))*S_in(i,j,1)
      enddo
      do k=2,nz-1 ; do i=is,ie
        ent(i,K+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
        c1(i,k) = ent(i,K) * b1(i)
        b1(i) = 1.0 / ((h(i,j,k) + d1(i)*ent(i,K)) + ent(i,K+1))
        d1(i) = b1(i) * (h(i,j,k) + d1(i)*ent(i,K))
        T_f(i,j,k) = b1(i) * (h(i,j,k)*T_in(i,j,k) + ent(i,K)*T_f(i,j,k-1))
        S_f(i,j,k) = b1(i) * (h(i,j,k)*S_in(i,j,k) + ent(i,K)*S_f(i,j,k-1))
      enddo ; enddo
      do i=is,ie
        c1(i,nz) = ent(i,nz) * b1(i)
        b1(i) = 1.0 / (h(i,j,nz) + d1(i)*ent(i,nz))
        T_f(i,j,nz) = b1(i) * (h(i,j,nz)*T_in(i,j,nz) + ent(i,nz)*T_f(i,j,nz-1))
        S_f(i,j,nz) = b1(i) * (h(i,j,nz)*S_in(i,j,nz) + ent(i,nz)*S_f(i,j,nz-1))
      enddo
      do k=nz-1,1,-1 ; do i=is,ie
        T_f(i,j,k) = T_f(i,j,k) + c1(i,k+1)*T_f(i,j,k+1)
        S_f(i,j,k) = S_f(i,j,k) + c1(i,k+1)*S_f(i,j,k+1)
      enddo ; enddo
    enddo
  endif

end subroutine vert_fill_TS


subroutine thickness_diffuse_init(Time, G, param_file, diag, CS)
  type(time_type),       intent(in) :: Time
  type(ocean_grid_type), intent(in) :: G
  type(param_file_type), intent(in) :: param_file
  type(diag_ptrs), target, intent(inout) :: diag
  type(thickness_diffuse_CS),     pointer    :: CS
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure containing pointers to common diagnostic fields.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
  integer :: nkbl

  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod = "GOLD_thickness_diffuse" ! This module's name.
  character(len=48)  :: flux_units

  if (associated(CS)) then
    call GOLD_error(WARNING, &
      "Thickness_diffuse_init called with an associated control structure.")
    return
  else ; allocate(CS) ; endif

  CS%diag => diag

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mod, version, tagname, "")
  call get_param(param_file, mod, "THICKNESSDIFFUSE", CS%thickness_diffuse, &
                 "If true, interfaces heights are diffused with a \n"//&
                 "coefficient of KHTH.", default=.false.)
  call get_param(param_file, mod, "KHTH", CS%Khth, &
                 "The background horizontal thickness diffusivity.", &
                 units = "m2 s-1", default=0.0)
  call get_param(param_file, mod, "KHTH_SLOPE_CFF", CS%KHTH_Slope_Cff, &
                 "The nondimensional coefficient in the Visbeck formula \n"//&
                 "for the interface depth diffusivity", units="nondim", &
                 default=0.0)
  call get_param(param_file, mod, "KHTH_MIN", CS%KHTH_Min, &
                 "The minimum horizontal thickness diffusivity.", &
                 units = "m2 s-1", default=0.0)
  call get_param(param_file, mod, "KHTH_MAX", CS%KHTH_Max, &
                 "The maximum horizontal thickness diffusivity.", &
                 units = "m2 s-1", default=0.0)
  call get_param(param_file, mod, "DETANGLE_INTERFACES", CS%detangle_interfaces, &
                 "If defined add 3-d structured enhanced interface height \n"//&
                 "diffusivities to horizonally smooth jagged layers.", &
                 default=.false.)
  CS%detangle_time = 0.0
  if (CS%detangle_interfaces) &
    call get_param(param_file, mod, "DETANGLE_TIMESCALE", CS%detangle_time, &
                 "A timescale over which maximally jagged grid-scale \n"//&
                 "thickness variations are suppressed.  This must be \n"//&
                 "longer than DT, or 0 to use DT.", units = "s", default=0.0) 
  call get_param(param_file, mod, "KHTH_SLOPE_MAX", CS%slope_max, &
                 "A slope beyond which the calculated isopycnal slope is \n"//&
                 "not reliable and is scaled away. This is used with \n"//&
                 "FULL_THICKNESSDIFFUSE.", units="nondim", default=0.01)
  call get_param(param_file, mod, "KD_SMOOTH", CS%kappa_smooth, &
                 "A diapycnal diffusivity that is used to interpolate \n"//&
                 "more sensible values of T & S into thin layers.", &
                 default=1.0e-6)
  call get_param(param_file, mod, "DEBUG", CS%debug, &
                 "If true, write out verbose debugging data.", default=.false.)


  if (G%Boussinesq) then ; flux_units = "meter3 second-1"
  else ; flux_units = "kilogram second-1" ; endif

  CS%id_uhGM = register_diag_field('ocean_model', 'uhGM', G%axesul, Time, &
           'Time Mean Diffusive Zonal Thickness Flux', flux_units)
  if (CS%id_uhGM > 0) call safe_alloc_ptr(diag%uhGM,G%Isdq,G%Iedq,G%jsd,G%jed,G%ke)
  CS%id_vhGM = register_diag_field('ocean_model', 'vhGM', G%axesvl, Time, &
           'Time Mean Diffusive Meridional Thickness Flux', flux_units)
  if (CS%id_vhGM > 0) call safe_alloc_ptr(diag%vhGM,G%isd,G%ied,G%Jsdq,G%Jedq,G%ke)
  CS%id_GMwork = register_diag_field('ocean_model', 'GMwork', G%axesh1, Time, &
           'Time Mean Integral Work done by Diffusive Thickness Flux', 'Watt meter-2')
  if (CS%id_GMwork > 0) call safe_alloc_ptr(CS%GMwork,G%isd,G%ied,G%jsd,G%jed)
  CS%id_KH_u = register_diag_field('ocean_model', 'KHTH_u', G%axesui, Time, &
           'Thickness Diffusivity at U-point', 'meter second-2')
  CS%id_KH_v = register_diag_field('ocean_model', 'KHTH_v', G%axesvi, Time, &
           'Thickness Diffusivity at V-point', 'meter second-2')
  CS%id_KH_u1 = register_diag_field('ocean_model', 'KHTH_u1', G%axesu1, Time, &
           'Thickness Diffusivity at U-points (2-D)', 'meter second-2')
  CS%id_KH_v1 = register_diag_field('ocean_model', 'KHTH_v1', G%axesv1, Time, &
           'Thickness Diffusivity at V-points (2-D)', 'meter second-2')

 ! CS%id_sfn_x =  register_diag_field('ocean_model', 'sfn_x', G%axesui, Time, &
 !          'Parameterized Zonal Overturning Streamfunction', 'meter3 second-1')
 ! CS%id_sfn_y =  register_diag_field('ocean_model', 'sfn_y', G%axesvi, Time, &
 !          'Parameterized Meridional Overturning Streamfunction', 'meter3 second-1')
 ! CS%id_sfn_slope_x =  register_diag_field('ocean_model', 'sfn_sl_x', G%axesui, Time, &
 !          'Parameterized Zonal Overturning Streamfunction from Interface Slopes', 'meter3 second-1')
 ! CS%id_sfn_slope_y =  register_diag_field('ocean_model', 'sfn_sl_y', G%axesvi, Time, &
 !          'Parameterized Meridional Overturning Streamfunction from Interface Slopes', 'meter3 second-1')

end subroutine thickness_diffuse_init

subroutine thickness_diffuse_end(CS)
  type(thickness_diffuse_CS), pointer :: CS
  if(associated(CS)) deallocate(CS)
end subroutine thickness_diffuse_end

end module GOLD_thickness_diffuse
